Preamble

This vignette is specifically focused on building references needed to stage samples with RAPToR. For a more general use of the package, see the "RAPToR" vignette.

Building references is one of the key aspects of RAPToR. An appropriate reference is needed to stage samples, and since they are re-usable we believe they are worth the trouble to set up.

Throughout this vignette, you will see the general workflow of building a reference from the selection of an appropriate dataset, to choosing and validating a model for interpolation. In the midst of the explanations, examples will be given using the same dsaeschimann2017 and dshendriks2014 datasets as in the general usage vignette ( you can check out how these were built at the end of the vignette ).

Finally, a few more examples of reference building on different datasets will be included at the end of this document.

I hope this material will be sufficient for your reference-building needs.

The data

Selecting / Preparing a dataset

Without a transcriptomic timecourse dataset spanning the developmental stages of your samples’ organism, I’m afraid there’s not much we can do for you ! Thankfully, these time-series experiments are (increasingly) numerous in the literature and most model organisms will have some kind of dataset we can use. You may also have your own time-series data on hand.

Databases

There are a few databases you can download datasets from. The most well-known are the Gene Expression Omnibus (GEO) and the Array Express.

Both of these databases have APIs to get their data from R (e.g the GEOquery package, as shown in the example dataset loading scripts).

What to look out for

Several points of the experimental design should be kept in mind when selecting a dataset for a reference.

Are there replicates ? If so, good. This means you can confirm the dynamics in your data are not noise. I would choose a sparser time-course with replicates over a higher-resolution experiment with one batch to build a reference.
Is the timepoint sampling even ? Transcriptomic profiling is expensive. Time courses experiments usually account for dynamic ranges of development (i.e, early or fast-changing stages are more sampled). For our purpose, we ideally want even sampling. However, if a dataset’s sampling respects dynamic ranges, you can still use it for reference building with a trick using ranks.
What’s the developmental range ? The bigger, the better ! (Though as long as the samples to stage fit within, it is enough).
Is the profiling done on whole-organism or specific parts ? You should aim for whole-organism profiling. Even if you have to stage tissue-specific samples, a good quality reference should still allow you to get accurate estimates, perhaps restraining the geneset used for staging.
What’s the profiling technology ? RNA-seq is much, much cleaner than MicroArray data. Sometimes, you just have to make do. When using RNA-seq counts, apply some within-sample normalization (e.g. TPM) to reflect gene expression accurately.

ID formatting

One of the plagues of bioinformatics is the large and fast-changing set of IDs for genes, transcripts, etc. When you build a reference, you should always convert it to IDs that are conventional and stable. We like to use the organism-specific IDs (e.g, Wormbase for C. elegans : WBGene00016153, Flybase for Drosophila : FBgn0010583).

The ensembl biomart or its associated R package biomaRt are a very useful ressource to get gene, transcript or probe ID lists for conversion.

Below is a code snippet with an example getting gene IDs for drosophila with biomaRt.

requireNamespace("biomaRt", quietly = TRUE)

# setup connection to ensembl
mart <- biomaRt::useMart("ensembl", dataset = "dmelanogaster_gene_ensembl")

# get list of attributes
droso_genes <- biomaRt::getBM(attributes = c("ensembl_gene_id", 
                                             "ensembl_transcript_id",
                                             "external_gene_name",
                                             "flybase_gene_id"),
                              mart = mart)

head(droso_genes)
#>   ensembl_gene_id ensembl_transcript_id external_gene_name flybase_gene_id
#> 1     FBgn0053882           FBtr0091886      His2B:CG33882     FBgn0053882
#> 2     FBgn0035648           FBtr0113148            CG13288     FBgn0035648
#> 3     FBgn0035648           FBtr0331563            CG13288     FBgn0035648
#> 4     FBgn0032726           FBtr0342927            CG10621     FBgn0032726
#> 5     FBgn0032726           FBtr0342926            CG10621     FBgn0032726
#> 6     FBgn0032726           FBtr0081191            CG10621     FBgn0032726

When multiple probe or transcript IDs match a single gene ID, we usually go for mean-aggregation of expression values. This is taken care of with the format_ids() function.

Normalize and log expression

It’s common practice to normalize expression datasets (e.g. to account for technical bias). You may deal with many different profiling technologies when building references, and may join multiple datasets together for a reference.

To stay as consistent as possible, we apply quantile-normalization on our datasets regardless of its source. For this, we use the normalizeBetweenArrays() function of the limma package.

We also apply a \(log(X + 1)\).

dsaeschimann2017$g <- limma::normalizeBetweenArrays(dsaeschimann2017$g, method = "quantile")
dsaeschimann2017$g <- log1p(dsaeschimann2017$g)

dshendriks2014$g <- limma::normalizeBetweenArrays(dshendriks2014$g, method = "quantile")
dshendriks2014$g <- log1p(dshendriks2014$g)

Observing the data

It’s usually a good practice to take a look at what’s inside your data before doing anything else to it.

dsaeschimann2017$g[1:5,1:4]
#>                let.7.n2853._18hr let.7.n2853._20hr let.7.n2853._22hr let.7.n2853._24hr
#> WBGene00007063          7.501850         10.988212          10.45480          7.994587
#> WBGene00007064          8.023767          8.655388          14.21012          9.759401
#> WBGene00007065         15.919452         16.875057          15.23932         18.847718
#> WBGene00003525          1.416181         10.938876          13.42202          2.488798
#> WBGene00007067          1.765342          1.775650           2.77224          2.200257

head(dsaeschimann2017$p, n = 5)
#>                        title geo_accession           organism_ch1       strain
#> GSM2113587 let.7.n2853._18hr    GSM2113587 Caenorhabditis elegans let-7(n2853)
#> GSM2113588 let.7.n2853._20hr    GSM2113588 Caenorhabditis elegans let-7(n2853)
#> GSM2113589 let.7.n2853._22hr    GSM2113589 Caenorhabditis elegans let-7(n2853)
#> GSM2113590 let.7.n2853._24hr    GSM2113590 Caenorhabditis elegans let-7(n2853)
#> GSM2113591 let.7.n2853._26hr    GSM2113591 Caenorhabditis elegans let-7(n2853)
#>            time in development:ch1 age
#> GSM2113587                18 hours  18
#> GSM2113588                20 hours  20
#> GSM2113589                22 hours  22
#> GSM2113590                24 hours  24
#> GSM2113591                26 hours  26

Correlation

With our time series data, we can look at correlation heatmaps or boxplots of the samples to catch potential ouliers and observe the clear correlation between samples of similar development.

cor_dsaeschimann2017 <- cor(dsaeschimann2017$g, method = "spearman")

Plots

Code

# Heatmap
ord <- order(dsaeschimann2017$p$age)
heatmap(cor_dsaeschimann2017[ord, ord], Colv = NA, Rowv = NA, scale = "none", keep.dendro = F, margins = c(1,1),
        RowSideColors = as.numeric(dsaeschimann2017$p$strain[ord]), labRow = "", labCol = "")

par(xpd = T) # text may have to be tweaked to plot size
mtext(text = unique(dsaeschimann2017$p$age), side = 1, line = 4, at = seq(-.1,1.05, l = 11)) 

# Boxplot
boxplot(cor_dsaeschimann2017 ~ interaction(dsaeschimann2017$p$strain, dsaeschimann2017$p$age), col = 1:4, xaxt = "n", 
        ylab = "Spearman correlation", xlab = "age", at = seq(1,44, l = 55)[c(T,T,T,T,F)])
axis(side = 1, at = seq(2,42, l = 11), labels = unique(dsaeschimann2017$p$age))

legend(23,.86, fill = 1:4, legend = c("let-7", "lin-41", "let-7/lin-41", "N2"), bty = "n")

Plotting components

You can also plot components (PCA or ICA) with respect to time, to get an idea of the dynamics. For PCA, we want to perform a non-scaled, centered PCA. Centering must be done gene-wise, not sample-wise (hence the matrix rotation below).

pca_dsaeschimann2017 <- stats::prcomp(t(dsaeschimann2017$g), rank = 25,
                                      center = TRUE, scale = FALSE)

Plot

Code

par(mfrow = c(2,4))
invisible(
sapply(seq_len(8), function(i){
  plot(dsaeschimann2017$p$age, pca_dsaeschimann2017$x[,i], lwd = 2, col = dsaeschimann2017$p$strain,
       xlab = "age", ylab = "PC", main = paste0("PC", i))
  
  # connect the dots
  sapply(seq_along(levels(dsaeschimann2017$p$strain)), function(l){
    s <- which(dsaeschimann2017$p$strain == levels(dsaeschimann2017$p$strain)[l])
    points(dsaeschimann2017$p$age[s], pca_dsaeschimann2017$x[s,i], col = l, 
           type = 'l', lty = 2)
  })

  if(i == 1)
    legend("topleft", bty = 'n', legend = c("let-7", "lin-41", "let-7/lin-41", "N2"),
           pch = c(rep(1, 4)), lty = c(rep(NA, 4)), col = c(1:4), lwd = 3)
})
)

In this data, we see that between different strains, we get very consistent dynamics. Also, PC2 and PC3 capture an oscillatory dynamic which is characteristic of the molting processes of C. elegans larval development.

Plotting random genes

Another approach can be to look at a few random genes. You get a first hand look at the noise in your data.

Plots

Code

set.seed(10)
gtp <- sample(nrow(dsaeschimann2017$g), size = 4)

par(mfrow = c(1,4))
invisible(
sapply(gtp, function(i){
  plot(dsaeschimann2017$p$age, dsaeschimann2017$g[i,], lwd = 2, col = dsaeschimann2017$p$strain,
       xlab = "age", ylab = "GExpr", main = rownames(dsaeschimann2017$g)[i])
  
  # connect the dots
  sapply(seq_along(levels(dsaeschimann2017$p$strain)), function(l){
    s <- which(dsaeschimann2017$p$strain == levels(dsaeschimann2017$p$strain)[l])
    points(dsaeschimann2017$p$age[s], dsaeschimann2017$g[i,s], col = l, 
           type = 'l', lty = 2)
  })

  if(i == gtp[4])
    legend("topleft", bty = 'n', legend = c("let-7", "lin-41", "let-7/lin-41", "N2"),
           pch = c(rep(1, 4)), lty = c(rep(NA, 4)), col = c(1:4), lwd = 3)
})
)

The gene expression interpolation model

To increase the resolution of our time series, we are faced with a very unbalanced regression problem. We essentially want to predict tens of thousands of dependent variables (genes) with our few independent variables (time, batch, …).

We refer to the gene expression interpolation model as GEIM in the following text.

About the model

The principal strategy we put forward for predicting on such a large scale of output variables, is to interpolate in a dimensionally reduced space. We propose to do this on Principal Components or Independant Components ( Independant Component Analysis ).

Both PCA and ICA perform the same type of linear transformation on the data, they just maximize a different criteria. PCA maximizes the variance of each component and ICA their independance. We get the following :

\[ X_{(m\times n)} = G_{(m\times c)}S^{T}_{(n\times c)} \] with \(X\), the matrix of \(m\) genes by \(n\) samples, \(G\) the gene loadings (\(m\) genes by \(c\) components) and \(S^T\) the sample loadings (\(n\) samples by \(c\) components). \(S\) is what’s usually looked at when performing a PCA (or ICA) on gene expression data, to look at the samples in the component space. It’s what we plotted in the section on observing data, for instance.

Alter, Brown, and Botstein (2000) previously demonstrated that singular value decomposition of gene expression data can be taken as “eigengenes”, giving a global picture of the dynamics of gene expression. We essentially use the same property for a GEIM. We build a model on the columns of \(S^T\) (eigengenes), predict in the component space, and reconstruct the gene expression data by a matrix product with the gene loadings.

We’ve implemented 2 model types : Generalized Additive Models (GAMs, the default) and Generalized Linear Models (GLMs). GAMs rely on the gam() function of the mgcv package, and GLMs on the glm() function of the stats core package.

As you’ll see in the next section, a standard R formula will be specifed to the model. This formula can make use of all the tools one can use with gam() or glm(), most notably the variety of polynomial or smoothing splines implemented through the s() function for GAMs. With GLMs, you can also use splines from the splines core package, such as ns() for natural cubic splines.

The GEIM interface

GEIMs, are built with the ge_im() function, which outputs a geim object. This function takes as input 3 key arguments :

X : the gene expression matrix of your time series (genes as rows, samples as columns)
p : a dataframe of pheno data, samples as rows. This should include the age/time variable and any other covariates you want to include in the model (e.g batch, strain)
formula : the model formula. This should be a standard R formula, using terms found in p. It should start with X ~.

Another important argument is the number of components used for the interpolation, nc.

For example, using the dsaeschimann2017 dataset we could build the following model.

m_dsaeschimann2017 <- ge_im(X = dsaeschimann2017$g, p = dsaeschimann2017$p, 
                            formula = "X ~ s(age, bs = 'ts') + strain", nc = 32)

Note that a single model formula is specified and applied to all the components, but the models are fitted independantly on each component.

Feel free to have a look at the documentation of the function for additionnal parameters ?ge_im().

To get model predictions, you simply use the predict() function, like for any standard R model.

Finding the appropriate model and parameters

Model type

There are 5 types of GEIM you can fit with the ge_im() function.

A GAM on PCA components (method = "gam", dim_red = "pca") (default)
A GLM on PCA components (method = "glm", dim_red = "pca")
A GAM on ICA components (method = "gam", dim_red = "ica")
A GLM on ICA components (method = "glm", dim_red = "ica")
A gene-by-gene linear model directly on the gene expression matrix (method = "limma")

The default option is a robust choice when applying a smoothing spline to the data.

Using PCA or ICA components yields near-identical results in most scenarios. ICA tends to outperform PCA when the data is very noisy. This is by design since ICA essentially performs signal extraction. It is however slower, especially if nc is large.

The last option ("limma") corresponds to a solution that makes no effort to reduce the dimensionality of the problem (the dim_red and nc arguments are ignored). As a result, there is no information loss or bias introduced by dimension reduction. This approach is however very sensitive to noise. The model is fit with the lmFit() function of the limma package (hence the name).

Note that when using GAMs, there is by definition no interaction possible between terms. It is possible to include a by argument to the s() function, but this essentialy corresponds to separate fits on the levels of the specified factor/variable.

Model performance

One can use a number of criteria to evaluate model performance. We provide the mperf() function to compute the indices described below by inputing the data and the model predictions.

In the formulas below, \(X\) corresponds to the input gene expression matrix (\(m\) genes as rows, \(n\) samples as columns), \(\hat{X}\) to the model predictions. \(x_i\) corresponds to row \(i\) of matrix \(X\) and \(x_i^{(j)}\) to sample \(j\) of that row. This notation is derived from the general regression problem, where \(X^T\) corresponds to the set of \(m\) dependant variable to predict.

aCC : average Correlation Coefficient.

\[ aCC = \frac{1}{m}\sum^{m}_{i=1}{CC} = \frac{1}{m}\sum^{m}_{i=1}{\cfrac{\sum^{n}_{j=1}{(x_i^{(j)}-\bar{x}_i)(\hat{x}_i^{(j)}-\bar{\hat{x}}_i)}}{\sqrt{\sum^{n}_{j=1}{(x_i^{(j)}-\bar{x}_i)^2(\hat{x}_i^{(j)}-\bar{\hat{x}}_i)^2}}}} \]

aRE : average Relative Error.

\[ a\delta = \frac{1}{m}\sum^{m}_{i=1}{\delta} = \frac{1}{m} \sum^{m}_{i=1} \frac{1}{n} \sum^{n}_{j=1} \cfrac{| x_i^{(j)} - \hat{x}_i^{(j)} | }{x_i^{(j)}} \]

MSE : Mean Squared Error.

\[ MSE = \frac{1}{m} \sum^{m}_{i=1} \frac{1}{n} \sum^{n}_{j=1} (x_i^{(j)} - \hat{x}_i^{(j)} )^2 \]

aRMSE : average Root MSE.

\[ aRMSE = \frac{1}{m}\sum^{m}_{i=1}{RMSE} = \frac{1}{m} \sum^{m}_{i=1} \sqrt{\cfrac{\sum^{n}_{j=1} (x_i^{(j)} - \hat{x}_i^{(j)} )^2}{n}} \]

Note that these indices are computed and averaged with respect to variables, not observations. You can either get the overall (averaged) index or the value per gene with the global parameter.

g_mp <- mperf(dsaeschimann2017$g, predict(m_dsaeschimann2017), is.t = T)
g_mp
#> $aCC
#> [1] 0.7977299
#> 
#> $aRE
#> [1] 0.1301014
#> 
#> $MSE
#> [1] 0.01431891
#> 
#> $aRMSE
#> [1] 0.1196617
ng_mp <- mperf(dsaeschimann2017$g, predict(m_dsaeschimann2017), is.t = T, global = F)
ng_mp <- lapply(ng_mp, na.omit) # remove NAs (eg. 0 variance genes)
ng_mp$aRE <- ng_mp$aRE[ng_mp$aRE < Inf] # remove Inf values (/0)

It’s possible to check the model performance by looking at the index distributions over all genes, e.g. :

Plots

Code

par(mfrow = c(2,2))
invisible(
sapply(names(ng_mp), function(idx){
  rg <- range(na.omit(ng_mp[[idx]]))
  
  # label position
  if(idx == "aCC"){
    pos <- 2
  } else {
    pos <- 4
  }
  # estimate density curve
  d <- density(na.omit(ng_mp[[idx]]), from = rg[1], to = rg[2])
  
  plot(d, main = paste0(gsub("a", "", idx, fixed = T), " density (", length(ng_mp[[idx]]), " genes)"), 
       xlab = idx, lwd = 2)
  # display global value
  abline(v = g_mp[[idx]], lty = 2, lwd = 2, col = "firebrick")
  text(g_mp[[idx]], .9*max(d$y), pos = pos, labels = idx, font = 2, col = "firebrick")
})
)

Number of components

The number of components to use for the interpolation is by default set to the number of samples. However, we recommend to set a cutoff on explained variance of PCA components to select it. You can do this (on PCA explained variance) even if you choose to use ICA for interpolation.

For example, on the dsaeschimann2017 dataset, we set the threshold at \(99\%\) :

nc <- sum(summary(pca_dsaeschimann2017)$importance[3,] < .99) + 1
nc
#> [1] 32

Note that this threshold must be set in accordance with the noise in the data. For example, in very noisy data, would you consider that \(99\%\) of the variance in the dataset corresponds to meaningful dynamics ?

You can also define nc by plotting your components and stopping after the components stop capturing meaningful variation (dynamics) with respect to time/age.

In a very noisy dataset, you may have to keep very few components (\(<5\)) for the interpolation.

Comparing formulas

Choosing from different splines (and/or parameters) can be done with cross-validation (CV) through the use of the ge_imCV() function. The function inputs the X, p and a formula_list to test. Other parameters on the CV itself can also be given (e.g. training set size).

The default training/validation set ratio is cv.s = 0.8, so \(80\%\) of the data is used to build the model. When including (factor) covariates in the model, the training set is built such that all groups are proportionately represented in the training set (based on terms of the first formula in the list). The number of repeats to do for the CV is defined by cv.n.

Note that the model type (GAM/GLM, PCA/ICA) is fixed for all formulas in one call of ge_imCV().

ge_imCV() computes the indices of model performance with mperf(), excluding aCC due to computing time. Indices are computed on the validation set (CV Error) and on the training set (Model PerFormance).

Below is an example of usage to choose between 4 available GAM smooth terms on the dsaeschimann2017 GEIM.

smooth_methods <- c("tp", "ts", "cr", "ps")
flist <- as.list(paste0("X ~ s(age, bs = \'", smooth_methods, "\') + strain"))
flist
#> [[1]]
#> [1] "X ~ s(age, bs = 'tp') + strain"
#> 
#> [[2]]
#> [1] "X ~ s(age, bs = 'ts') + strain"
#> 
#> [[3]]
#> [1] "X ~ s(age, bs = 'cr') + strain"
#> 
#> [[4]]
#> [1] "X ~ s(age, bs = 'ps') + strain"

cv_dsaeschimann2017 <- ge_imCV(X = dsaeschimann2017$g, p = dsaeschimann2017$p, formula_list = flist,
                  cv.n = 20, nc = nc, nb.cores = 3)
#> CV on 4 models. cv.n = 20 | cv.s = 0.8
#> 
#> ...Building training sets
#> ...Setting up cluster
#> ...Running CV
#> ...Cleanup and formatting

plot(cv_dsaeschimann2017, names = paste0("bs = ", smooth_methods), outline = F,
     swarmargs = list(cex = .8))

From the plots above, we can see the different splines perform similarly in this scenario, all could work. We chose ts (a thin-plate regression spline), as it seems to minimize CV error without much impact model performance.

You can also specify extra spline parameters for a fit. With s(), you can give the spline’s basis dimension (~= number of knots) to use with the k parameter. By default, the spline is a penalized spline, so it will not necessarily use k knots, but it will stay below that value. In our experience, this is not really necessary as the parameter estimation done with gam() is usually sufficient.

To force a spline with k basis dimension, you must also set fx = TRUE (note this fits a spline of dimension k on all components, whereas the penalized spline will adjust). You can look at the s() or choose.k documentation for further information.

Below is an example with the dsaeschimann2017 data.

ks <- c(4,6,8,10)
flistk <- as.list(c(
  "X ~ s(age, bs =  'cr') + strain",
  paste0("X ~ s(age, bs =  'cr', k = ", ks , ") + strain"), 
  paste0("X ~ s(age, bs =  'cr', k = ", ks , ", fx=TRUE) + strain")
  ))
flistk
#> [[1]]
#> [1] "X ~ s(age, bs =  'cr') + strain"
#> 
#> [[2]]
#> [1] "X ~ s(age, bs =  'cr', k = 4) + strain"
#> 
#> [[3]]
#> [1] "X ~ s(age, bs =  'cr', k = 6) + strain"
#> 
#> [[4]]
#> [1] "X ~ s(age, bs =  'cr', k = 8) + strain"
#> 
#> [[5]]
#> [1] "X ~ s(age, bs =  'cr', k = 10) + strain"
#> 
#> [[6]]
#> [1] "X ~ s(age, bs =  'cr', k = 4, fx=TRUE) + strain"
#> 
#> [[7]]
#> [1] "X ~ s(age, bs =  'cr', k = 6, fx=TRUE) + strain"
#> 
#> [[8]]
#> [1] "X ~ s(age, bs =  'cr', k = 8, fx=TRUE) + strain"
#> 
#> [[9]]
#> [1] "X ~ s(age, bs =  'cr', k = 10, fx=TRUE) + strain"

cv_dsaeschimann2017k <- ge_imCV(X = dsaeschimann2017$g, p = dsaeschimann2017$p, formula_list = flistk,
                   cv.n = 20, nc = nc, nb.cores = 3)
#> CV on 9 models. cv.n = 20 | cv.s = 0.8
#> 
#> ...Building training sets
#> ...Setting up cluster
#> ...Running CV
#> ...Cleanup and formatting

Predicting new timepoints

The entire reason we need GEIMs is to interpolate. We want a higher precision time series to stage samples.

To do this, as for any R model, we will use the predict() function. We first need to set up the new data to predict from.

n.inter <- 100 # nb of new timepoints
newdat <- data.frame(
  age = seq(min(dsaeschimann2017$p$age), max(dsaeschimann2017$p$age), l = n.inter),
  strain = rep("N2", n.inter) # we want to predict as N2 
  )
head(newdat)
#>        age strain
#> 1 18.00000     N2
#> 2 18.20202     N2
#> 3 18.40404     N2
#> 4 18.60606     N2
#> 5 18.80808     N2
#> 6 19.01010     N2

# predict 
pred_m_dsaeschimann2017 <- predict(m_dsaeschimann2017, newdata = newdat)

You can choose to output the predictions directly as the full gene expression matrix (default) or as components. Predictions as components with the as.c = TRUE option can be useful when checking the interpolation (see the next section).

pred_m_dsaeschimann2017_comp <- predict(m_dsaeschimann2017, newdata = newdat, as.c = TRUE)

Validating / Checking model predictions

We have built our model and predicted new data, now we can have a look at the interpolation results. To do this, one can do two things.

Checking model predictions against components

Checking model predictions against components allows you to see immediately if some dynamics get mishandled by the model.

Don’t hope for perfect fits ! It’s perfectly acceptable to have slight offsets. You may also notice some noisy components get “flattened”, with a null model fitted. These have little impact on the results as they usually correspond to a minuscule part of total variance. They also sometimes actually get rid of some unwanted variation.

On the dsaeschimann2017 dataset, the predictions of the first few components are plotted below.

Plots

Code

par(mfrow = c(2,4))
invisible(
sapply(seq_len(8), function(i){
  plot(dsaeschimann2017$p$age, pca_dsaeschimann2017$x[,i], lwd = 2, col = dsaeschimann2017$p$strain,
       xlab = "age", ylab = "PC", main = paste0("PC", i))
  
  # connect the dots
  sapply(seq_along(levels(dsaeschimann2017$p$strain)), function(l){
    s <- which(dsaeschimann2017$p$strain == levels(dsaeschimann2017$p$strain)[l])
    points(dsaeschimann2017$p$age[s], pca_dsaeschimann2017$x[s,i], col = l, 
           type = 'l', lty = 2)
  })
  
  # add model prediction
  points(newdat$age, pred_m_dsaeschimann2017_comp[, i], col = "royalblue", type = 'l', lwd = 3)
  
  if(i == 1)
    legend("topleft", bty = 'n', legend = c("let-7", "lin-41", "let-7/lin-41", "N2", "pred"),
           pch = c(rep(1, 4), NA), lty = c(rep(NA, 4), 1), col = c(1:4, "royalblue"), lwd = c(rep(3,4),4))
})
)

The interpolation should translates well back on the full expression matrix :

Plots

Code

par(mfrow = c(1,4))
invisible(
sapply(gtp, function(i){ # gtp is from the earlier gene plots
  plot(dsaeschimann2017$p$age, dsaeschimann2017$g[i,], lwd = 2, col = dsaeschimann2017$p$strain,
       xlab = "age", ylab = "GExpr", main = rownames(dsaeschimann2017$g)[i])
  
  # connect the dots
  sapply(seq_along(levels(dsaeschimann2017$p$strain)), function(l){
    s <- which(dsaeschimann2017$p$strain == levels(dsaeschimann2017$p$strain)[l])
    points(dsaeschimann2017$p$age[s], dsaeschimann2017$g[i,s], col = l, 
           type = 'l', lty = 2)
  })
  
  # add model prediction
  points(newdat$age, pred_m_dsaeschimann2017[i, ], col = "royalblue", type = 'l', lwd = 3)
  
  if(i == gtp[4])
    legend("topleft", bty = 'n', legend = c("let-7", "lin-41", "let-7/lin-41", "N2"),
           pch = c(rep(1, 4)), lty = c(rep(NA, 4)), col = c(1:4), lwd = 3)
})
)

Staging samples

Testing your reference on its final purpose is a good practice.

A first test is to stage the samples you used to build the reference on their own interpolated reference.

Finally, the best validation is to stage another time series from the literature on your reference if such data exists. This provides external validation, confirming the dynamics you interpolated on indeed correspond to development processes and not a specific effect of the dataset (which is unlikely, but not impossible).

We can do this with our example, and the dshendriks2014 dataset can be used for external validation.

# make a 'reference object' 
r_dsaeschimann2017 <- list(interpGE = pred_m_dsaeschimann2017, time.series = newdat$age)

ae_test_dsaeschimann2017 <- ae(dsaeschimann2017$g, r_dsaeschimann2017$interpGE, r_dsaeschimann2017$time.series)
ae_test_dshendriks2014 <- ae(dshendriks2014$g, r_dsaeschimann2017$interpGE, r_dsaeschimann2017$time.series)

Plots

Code

par(mfrow = c(1,2))
rg <- range(c(ae_test_dsaeschimann2017$age.estimates[,1], dsaeschimann2017$p$age))

# Plot 1
plot(ae_test_dsaeschimann2017$age.estimates[,1]~dsaeschimann2017$p$age, 
     xlab = "Chronological age", ylab = "Estimated age (dsaeschimann2017)", 
     xlim = rg, ylim = rg,
     main = "Chron. vs Estimated ages for dsaeschimann2017\n(on dsaeschimann2017 reference)", lwd = 2, 
     col = factor(dsaeschimann2017$p$strain))
# connect the dots
invisible(sapply(levels(factor(dsaeschimann2017$p$strain)), function(l){
  s <- dsaeschimann2017$p$strain == l
  points(ae_test_dsaeschimann2017$age.estimates[s,1]~dsaeschimann2017$p$age[s], type = 'l', 
         lty = 2, col = which(l==levels(factor(dsaeschimann2017$p$strain))))
}))
abline(a = 0, b = 1, lty = 3, lwd = 2) # x = y

legend("bottomright", legend = c("let-7", "lin-41", "let-7/lin-41", "N2", "x = y"), 
       lwd=3, col=c(1:4, 1), bty='n', pch = c(1,1,1,1,NA), lty = c(rep(NA, 4), 3))


# Plot 2
rg <- range(c(ae_test_dshendriks2014$age.estimates[,1], dshendriks2014$p$age))
plot(ae_test_dshendriks2014$age.estimates[,1]~dshendriks2014$p$age, 
     xlab = "Chronological age", ylab = "Estimated age (dsaeschimann2017)", 
     xlim = rg, ylim = rg,
     main = "Chron. vs Estimated ages for dshendriks2014\n(on dsaeschimann2017 reference)", lwd = 2)
# connect the dots
points(ae_test_dshendriks2014$age.estimates[,1] ~ dshendriks2014$p$age, type = 'l', lty = 2)
abline(a = 0, b = 1, lty = 3, lwd = 2) # x = y

legend("bottomright", legend = "x = y", lwd=3, col=1, lty = 3, bty='n')

Reference-Building examples

Here you will find a few examples of reference building on different organisms.

C. elegans larval development (used in the vignette above)
D. melanogaster embryonic development
Danio rerio embryonic development

C. elegans larval development

The data

The datasets in this example are those used for all the in-text examples throughout the reference-building vignette.

We are using the two C. elegans time series datasets.

A time series of larval development in 4 different strains published by Aeschimann et al. (2017), called dsaeschimann2017. This is the dataset used to build the reference. (Accession : GSE80157)
A high-resolution time series of late larval development published by Hendriks et al. (2014), called dshendriks2014. This is the dataset used for external validation. (Accession : GSE52861)

Code to generate dsaeschimann2017 and dshendriks2014 :

Hide

Show

Note : set the data_folder variable to an existing path on your system where you want to store the objects.

data_folder <- "../inst/extdata/"

requireNamespace("wormRef", quietly = T)
requireNamespace("utils", quietly = T)
requireNamespace("GEOquery", quietly = T) # May need to be installed with bioconductor
requireNamespace("Biobase", quietly = T)

raw2tpm <- function(rawcounts, genelengths){
  if(nrow(rawcounts) != length(genelengths))
    stop("genelengths must match nrow(rawcounts).")
  x <- rawcounts/genelengths
  return(t( t(x) * 1e6 / colSums(x) ))
}

fpkm2tpm <- function(fpkm){
  return(exp(log(fpkm) - log(colSums(fpkm)) + log(1e6)))
}

`dsaeschimann2017`

geo_dsaeschimann2017 <- "GSE80157"

g_url_dsaeschimann2017 <- GEOquery::getGEOSuppFiles(geo_dsaeschimann2017, makeDirectory = FALSE, fetch_files = FALSE)
g_file_dsaeschimann2017 <- paste0(data_folder, "dsaeschimann2017.txt.gz")
utils::download.file(url = as.character(g_url_dsaeschimann2017$url[2]), destfile = g_file_dsaeschimann2017)

X_dsaeschimann2017 <- read.table(gzfile(g_file_dsaeschimann2017), h=T, sep = '\t', stringsAsFactors = F, row.names = 1)

# convert to tpm & wb_id
X_dsaeschimann2017 <- X_dsaeschimann2017[rownames(X_dsaeschimann2017)%in%wormRef::Cel_genes$wb_id,]
X_dsaeschimann2017 <- raw2tpm(rawcounts = X_dsaeschimann2017, 
                              genelengths = wormRef::Cel_genes$transcript_length[match(rownames(X_dsaeschimann2017),
                                                                                       wormRef::Cel_genes$wb_id)])

# pheno data
P_dsaeschimann2017 <- Biobase::pData(GEOquery::getGEO(geo_dsaeschimann2017, getGPL = F)[[1]])
P_dsaeschimann2017[,10:34] <- NULL
P_dsaeschimann2017[, 3:8] <- NULL

colnames(P_dsaeschimann2017)[4] <- "strain"
P_dsaeschimann2017$strain <- factor(P_dsaeschimann2017$strain)
P_dsaeschimann2017$title <- make.names(P_dsaeschimann2017$title)

colnames(X_dsaeschimann2017) <- gsub('RNASeq_riboM_', '', colnames(X_dsaeschimann2017), fixed = T)
P_dsaeschimann2017$title <- gsub('RNASeq_riboM_', '', P_dsaeschimann2017$title, fixed = T)

# get age 
P_dsaeschimann2017$age <- as.numeric(sub('(\\d+)\\shours', '\\1', P_dsaeschimann2017$`time in development:ch1`))


X_dsaeschimann2017 <- X_dsaeschimann2017[, P_dsaeschimann2017$title]

dsaeschimann2017 <- list(g = X_dsaeschimann2017, p = P_dsaeschimann2017)
save(dsaeschimann2017, file = paste0(data_folder, "dsaeschimann2017.RData"), compress = "xz")

# cleanup
file.remove(g_file_dsaeschimann2017)
rm(geo_dsaeschimann2017, g_url_dsaeschimann2017, g_file_dsaeschimann2017, X_dsaeschimann2017, P_dsaeschimann2017)

`dshendriks2014`

geo_dshendriks2014 <- "GSE52861"

g_url_dshendriks2014 <- GEOquery::getGEOSuppFiles(geo_dshendriks2014, makeDirectory = FALSE, fetch_files = FALSE)
g_file_dshendriks2014 <- paste0(data_folder, "dshendriks2014.txt.gz")
utils::download.file(url = as.character(g_url_dshendriks2014$url[2]), destfile = g_file_dshendriks2014)

X_dshendriks2014 <- read.table(gzfile(g_file_dshendriks2014), h=T, sep = '\t', stringsAsFactors = F, row.names = 1)

# convert to tpm & wb_id
X_dshendriks2014 <- X_dshendriks2014[rownames(X_dshendriks2014)%in%wormRef::Cel_genes$wb_id,]
X_dshendriks2014 <- raw2tpm(rawcounts = X_dshendriks2014, 
                            genelengths = wormRef::Cel_genes$transcript_length[match(rownames(X_dshendriks2014),
                                                                                     wormRef::Cel_genes$wb_id)])


# pheno data
P_dshendriks2014 <- Biobase::pData(GEOquery::getGEO(geo_dshendriks2014, getGPL = F)[[1]])

# filter relevant fields/samples
P_dshendriks2014 <- P_dshendriks2014[(P_dshendriks2014$`strain:ch1` == 'N2') & (P_dshendriks2014$`growth protocol:ch1` == 'Continuous'), ]
P_dshendriks2014 <- P_dshendriks2014[, c("title", "geo_accession", "time in development:ch1")]

# get age 
P_dshendriks2014$age <- as.numeric(sub('(\\d+)\\shours', '\\1', P_dshendriks2014$`time in development:ch1`))


# formatting
P_dshendriks2014$title <- gsub('RNASeq_polyA_', '', 
                               gsub('hr', 'h', 
                                    gsub('-', '.', fixed = T, as.character(P_dshendriks2014$title))))
colnames(X_dshendriks2014) <- gsub('RNASeq_polyA_','', colnames(X_dshendriks2014))
X_dshendriks2014 <- X_dshendriks2014[, P_dshendriks2014$title]

dshendriks2014 <- list(g = X_dshendriks2014, p = P_dshendriks2014)
save(dshendriks2014, file = paste0(data_folder, "dshendriks2014.RData"), compress = "xz")

# cleanup
file.remove(g_file_dshendriks2014)
rm(geo_dshendriks2014, g_url_dshendriks2014, g_file_dshendriks2014, X_dshendriks2014, P_dshendriks2014)

Normalization & Quick look

dsaeschimann2017$g <- limma::normalizeBetweenArrays(dsaeschimann2017$g, method = "quantile")
dsaeschimann2017$g <- log1p(dsaeschimann2017$g)

dshendriks2014$g <- limma::normalizeBetweenArrays(dshendriks2014$g, method = "quantile")
dshendriks2014$g <- log1p(dshendriks2014$g)

dsaeschimann2017$g[1:5,1:4]
#>                let.7.n2853._18hr let.7.n2853._20hr let.7.n2853._22hr let.7.n2853._24hr
#> WBGene00007063         2.1206604          2.469532          2.373273          2.175924
#> WBGene00007064         2.1621558          2.260804          2.661102          2.354485
#> WBGene00007065         2.7763061          2.847833          2.727037          2.960098
#> WBGene00003525         0.9434159          2.466223          2.609585          1.313603
#> WBGene00007067         1.0787531          1.081964          1.350796          1.236899

head(dsaeschimann2017$p, n = 5)
#>                        title geo_accession           organism_ch1       strain
#> GSM2113587 let.7.n2853._18hr    GSM2113587 Caenorhabditis elegans let-7(n2853)
#> GSM2113588 let.7.n2853._20hr    GSM2113588 Caenorhabditis elegans let-7(n2853)
#> GSM2113589 let.7.n2853._22hr    GSM2113589 Caenorhabditis elegans let-7(n2853)
#> GSM2113590 let.7.n2853._24hr    GSM2113590 Caenorhabditis elegans let-7(n2853)
#> GSM2113591 let.7.n2853._26hr    GSM2113591 Caenorhabditis elegans let-7(n2853)
#>            time in development:ch1 age
#> GSM2113587                18 hours  18
#> GSM2113588                20 hours  20
#> GSM2113589                22 hours  22
#> GSM2113590                24 hours  24
#> GSM2113591                26 hours  26

Correlation Matrix

Plotting components

pca_dsaeschimann2017 <- stats::prcomp(t(dsaeschimann2017$g), rank = 25,
                                      center = TRUE, scale = FALSE)

Model fitting

Component number

nc <- sum(summary(pca_dsaeschimann2017)$importance[3,] < .99) + 1
nc
#> [1] 32

Model

m_dsaeschimann2017 <- ge_im(X = dsaeschimann2017$g, p = dsaeschimann2017$p, 
                            formula = "X ~ s(age, bs = 'ts') + strain", nc = nc)

#>       aCC       aRE        MSE     aRMSE
#> 0.7977299 0.1301014 0.01431891 0.1196617

Validation

Predict

# setup newdat
n.inter <- 100 # nb of new timepoints
newdat <- data.frame(
  age = seq(min(dsaeschimann2017$p$age), max(dsaeschimann2017$p$age), l = n.inter),
  strain = rep("N2", n.inter) # we want to predict as N2 
  )
head(newdat)
#>        age strain
#> 1 18.00000     N2
#> 2 18.20202     N2
#> 3 18.40404     N2
#> 4 18.60606     N2
#> 5 18.80808     N2
#> 6 19.01010     N2

# predict 
pred_m_dsaeschimann2017 <- predict(m_dsaeschimann2017, newdata = newdat)
pred_m_dsaeschimann2017_comp <- predict(m_dsaeschimann2017, newdata = newdat, as.c = TRUE)

Plot component predictions

Build reference & stage samples

# make a 'reference object' 
r_dsaeschimann2017 <- list(interpGE = pred_m_dsaeschimann2017, time.series = newdat$age)

ae_dsaeschimann2017 <- ae(dsaeschimann2017$g, r_dsaeschimann2017$interpGE, r_dsaeschimann2017$time.series)
ae_dshendriks2014 <- ae(dshendriks2014$g, r_dsaeschimann2017$interpGE, r_dsaeschimann2017$time.series)

D. melanogaster embryonic development

The data

We are using two Drosophila melanogaster embryonic development time series datasets. The dataset used to build the reference was chosen with a very low time resolution on purpose to display the effectiveness of interpolating on gene expression data

A time series Drosophila embryonic development, part of the modENCODE project and published by Graveley et al. (2011), hereafter called dsgraveley2011. This is the dataset used to build the reference. (Data downloaded from fruitfly.org)
A high-resolution time series of embryonic development published by Levin et al. (2016), called dslevin2016dmel. This is the dataset used for external validation. (Accession : GSE60471)

Code to generate dsgraveley2011 and dslevin2016dmel :

Hide

Show

Note : set the data_folder variable to an existing path on your system where you want to store the objects.

data_folder <- "../inst/extdata/"

requireNamespace("utils", quietly = T)
requireNamespace("GEOquery", quietly = T) # May need to be installed with bioconductor
requireNamespace("Biobase", quietly = T)

raw2tpm <- function(rawcounts, genelengths){
  if(nrow(rawcounts) != length(genelengths))
    stop("genelengths must match nrow(rawcounts).")
  x <- rawcounts/genelengths
  return(t( t(x) * 1e6 / colSums(x) ))
}

fpkm2tpm <- function(fpkm){
  return(exp(log(fpkm) - log(colSums(fpkm)) + log(1e6)))
}

requireNamespace("biomaRt", quietly = TRUE)

mart <- biomaRt::useMart("ensembl", dataset = "dmelanogaster_gene_ensembl")
droso_genes <- biomaRt::getBM(attributes = c("ensembl_gene_id", 
                                             "ensembl_transcript_id",
                                             "external_gene_name",
                                             "transcript_length"),
                              mart = mart)
colnames(droso_genes)[1:3] <- c("fb_id", "transcript_id", "gene_name")

rm(mart)

`dsgraveley2011`

g_url_dsgraveley2011 <- "ftp://ftp.fruitfly.org/pub/download/modencode_expression_scores/Celniker_Drosophila_Annotation_20120616_1428_allsamps_MEAN_gene_expression.csv.gz"
g_file_dsgraveley2011 <- paste0(data_folder, "dsgraveley2011.csv.gz")
utils::download.file(g_url_dsgraveley2011, destfile = g_file_dsgraveley2011)


X_dsgraveley2011 <- read.table(gzfile(g_file_dsgraveley2011), sep = ',', row.names = 1, h = T)

# convert gene ids to FBgn
X_dsgraveley2011 <- RAPToR::format_ids(X_dsgraveley2011, droso_genes, from = "gene_name", to = "fb_id")

# select embryo time series samples
X_dsgraveley2011 <- X_dsgraveley2011[,1:12]

P_dsgraveley2011 <- data.frame(sname = colnames(X_dsgraveley2011),
                               age = as.numeric(gsub("em(\\d+)\\.\\d+hr", "\\1", colnames(X_dsgraveley2011))),
                               stringsAsFactors = FALSE)

dsgraveley2011 <- list(g = X_dsgraveley2011, p = P_dsgraveley2011)

save(dsgraveley2011, file = paste0(data_folder, "dsgraveley2011.RData"), compress = "xz")

# cleanup
file.remove(g_file_dsgraveley2011)
rm(g_url_dsgraveley2011, g_file_dsgraveley2011, X_dsgraveley2011, P_dsgraveley2011)

`dslevin2016dmel`

geo_dslevin2016dmel <- "GSE60471"

g_url_dslevin2016dmel <- GEOquery::getGEOSuppFiles(geo_dslevin2016dmel, makeDirectory = FALSE, fetch_files = FALSE)
g_file_dslevin2016dmel <- paste0(data_folder, "dslevin2016dmel.txt.gz")
utils::download.file(url = as.character(g_url_dslevin2016dmel$url[3]), destfile = g_file_dslevin2016dmel)

X_dslevin2016dmel <- read.table(gzfile(g_file_dslevin2016dmel), h = T, sep = '\t', as.is = T, row.names = 1, comment.char = "")

# filter poor quality samples
cm_dslevin2016dmel <- RAPToR::cor.gene_expr(X_dslevin2016dmel, X_dslevin2016dmel)
f_dslevin2016dmel <- which(0.6 > apply(cm_dslevin2016dmel, 1, quantile, probs = .99))
X_dslevin2016dmel <- X_dslevin2016dmel[, -f_dslevin2016dmel]

# convert to tpm & FBgn

X_dslevin2016dmel <- X_dslevin2016dmel[rownames(X_dslevin2016dmel)%in%droso_genes$fb_id,]
X_dslevin2016dmel <- raw2tpm(rawcounts = X_dslevin2016dmel, 
                            genelengths = droso_genes$transcript_length[match(rownames(X_dslevin2016dmel),
                                                                              droso_genes$fb_id)])

# pheno data
P_dslevin2016dmel <- Biobase::pData(GEOquery::getGEO(geo_dslevin2016dmel, getGPL = F)[[1]])

# filter relevant fields/samples
P_dslevin2016dmel <- P_dslevin2016dmel[, c("title", "geo_accession", "time (minutes cellularization stage):ch1")]
colnames(P_dslevin2016dmel)[3] <- "time"
P_dslevin2016dmel$title <- as.character(P_dslevin2016dmel$title)

P_dslevin2016dmel <- P_dslevin2016dmel[P_dslevin2016dmel$title %in% colnames(X_dslevin2016dmel),]
X_dslevin2016dmel <- X_dslevin2016dmel[, P_dslevin2016dmel$title]

# formatting
P_dslevin2016dmel$title <- gsub('Metazome_Drosophila_timecourse_', '', P_dslevin2016dmel$title)
colnames(X_dslevin2016dmel) <- P_dslevin2016dmel$title

P_dslevin2016dmel$age <- as.numeric(P_dslevin2016dmel$time) / 60

dslevin2016dmel <- list(g = X_dslevin2016dmel, p = P_dslevin2016dmel)
save(dslevin2016dmel, file = paste0(data_folder, "dslevin2016dmel.RData"), compress = "xz")

# cleanup
file.remove(g_file_dslevin2016dmel)
rm(geo_dslevin2016dmel, g_url_dslevin2016dmel, g_file_dslevin2016dmel, 
   X_dslevin2016dmel, P_dslevin2016dmel, 
   cm_dslevin2016dmel, f_dslevin2016dmel)

rm(droso_genes, raw2tpm, fpkm2tpm)

Normalization & Quick look

dsgraveley2011$g <- limma::normalizeBetweenArrays(dsgraveley2011$g, method = "quantile")
dsgraveley2011$g <- log1p(dsgraveley2011$g)

dslevin2016dmel$g <- limma::normalizeBetweenArrays(dslevin2016dmel$g, method = "quantile")
dslevin2016dmel$g <- log1p(dslevin2016dmel$g)

dsgraveley2011$g[1:5, 1:5]
#>               em0.2hr   em2.4hr   em4.6hr   em6.8hr  em8.10hr
#> FBgn0000003 3.9651391 4.2738527 3.3174101 4.5644242 4.6982706
#> FBgn0000008 1.2949845 0.9215699 0.6958672 0.6476801 0.7445991
#> FBgn0000014 0.5099295 0.9512866 1.3952815 1.8610406 1.8421960
#> FBgn0000015 0.2435639 0.6423988 1.0511912 1.1094674 1.0194280
#> FBgn0000017 1.7968429 2.0901351 1.3389420 1.5336183 1.6777064


head(dsgraveley2011$p, n = 5)
#>      sname age
#> 1  em0.2hr   0
#> 2  em2.4hr   2
#> 3  em4.6hr   4
#> 4  em6.8hr   6
#> 5 em8.10hr   8

Correlation Matrix

Plotting components

pca_dsgraveley2011 <- stats::prcomp(t(dsgraveley2011$g), rank = 12,
                                    center = TRUE, scale = FALSE)

Model fitting

Component number

nc <- sum(summary(pca_dsgraveley2011)$importance[3,] < .99) + 1
nc
#> [1] 8

Model

m_dsgraveley2011 <- ge_im(X = dsgraveley2011$g, p = dsgraveley2011$p, formula = "X ~ s(age, bs = 'cr')", nc = nc)

#>        aCC       aRE         MSE      aRMSE
#>  0.9400184  0.389298 0.009916307 0.09958065

Validation

Predict

# setup newdat
n.inter <- 100 # nb of new timepoints
newdat <- data.frame(
  age = seq(min(dsgraveley2011$p$age), max(dsgraveley2011$p$age), l = n.inter)
  )
head(newdat)
#>         age
#> 1 0.0000000
#> 2 0.2222222
#> 3 0.4444444
#> 4 0.6666667
#> 5 0.8888889
#> 6 1.1111111


# predict 
pred_m_dsgraveley2011 <- predict(m_dsgraveley2011, newdata = newdat)
pred_m_dsgraveley2011_comp <- predict(m_dsgraveley2011, newdata = newdat, as.c = TRUE)

Plot component predictions

Build reference & stage samples

# make a 'reference object' 
r_dsgraveley2011 <- list(interpGE = pred_m_dsgraveley2011, time.series = newdat$age)

ae_dsgraveley2011 <- ae(dsgraveley2011$g, r_dsgraveley2011$interpGE, r_dsgraveley2011$time.series)
ae_dslevin2016dmel <- ae(dslevin2016dmel$g, r_dsgraveley2011$interpGE, r_dsgraveley2011$time.series)

Notice here, that our validation dataset’s estimates appear quite noisy. However, if we look at the dynamics of the dslevin2016dmel data, we’ll see that the chronological age specified for the samples is erroneous.

pca_dslevin2016dmel <- stats::prcomp(t(dslevin2016dmel$g), rank = 20,
                                     center = TRUE, scale = FALSE)

This demonstrates the difficulty of producing high-resolution time series due to developmental asynchronicity between the samples.

Danio rerio embryonic development

The data

We are using two Danio rerio (zebrafish) embryonic development time series datasets. The dataset used to build the reference has uneven time sampling, as can often be the case.

We show a trick using ranks to build an adequate model in order to avoid interpolation bias.

The datasets are

A time series of zebrafish embryonic development (with uneven sampling) published by White et al. (2017), hereafter called dswhite2017. This is the dataset used to build the reference. (Data accessible in the publication)
A high-resolution time series of embryonic development published by Levin et al. (2016), hereafter called dslevin2016zeb. This is the dataset used for external validation. (Accession : GSE60619)

Code to generate dswhite2017 and dslevin2016zeb :

Hide

Show

Note : set the data_folder variable to an existing path on your system where you want to store the objects.

data_folder <- "../inst/extdata/"

requireNamespace("utils", quietly = T)
requireNamespace("GEOquery", quietly = T) # May need to be installed with bioconductor
requireNamespace("Biobase", quietly = T)

raw2tpm <- function(rawcounts, genelengths){
  if(nrow(rawcounts) != length(genelengths))
    stop("genelengths must match nrow(rawcounts).")
  x <- rawcounts/genelengths
  return(t( t(x) * 1e6 / colSums(x) ))
}

fpkm2tpm <- function(fpkm){
  return(exp(log(fpkm) - log(colSums(fpkm)) + log(1e6)))
}

requireNamespace("biomaRt", quietly = TRUE)

mart <- biomaRt::useMart("ensembl", dataset = "drerio_gene_ensembl")
zeb_genes <- biomaRt::getBM(attributes = c("ensembl_gene_id", "transcript_length"), 
                            mart = mart)
rm(mart)

`dswhite2017`

p_url_dswhite2017 <- "http://europepmc.org/articles/PMC5690287/bin/elife-30860-supp1.tsv"

g_url_dswhite2017 <- "http://europepmc.org/articles/PMC5690287/bin/elife-30860-supp2.tsv"
g_file_dswhite2017 <- paste0(data_folder, "dswhite2017.tsv")
utils::download.file(g_url_dswhite2017, destfile = g_file_dswhite2017)

X_dswhite2017 <- read.table(g_file_dswhite2017, h = T, sep  ="\t", as.is = T, quote = "\"")
rownames(X_dswhite2017) <- X_dswhite2017$Gene.ID
X_dswhite2017 <- X_dswhite2017[,-(1:8)]


# convert to tpm & ensembl_id
X_dswhite2017 <- X_dswhite2017[rownames(X_dswhite2017)%in%zeb_genes$ensembl_gene_id,]
X_dswhite2017 <- raw2tpm(rawcounts = X_dswhite2017, 
                            genelengths = zeb_genes$transcript_length[match(rownames(X_dswhite2017),
                                                                            zeb_genes$ensembl_gene_id)])

# pheno data
P_dswhite2017 <- read.table(p_url_dswhite2017, h = T, sep = "\t", as.is = T)
P_dswhite2017 <- P_dswhite2017[P_dswhite2017$sequencing == "RNASeq", c("sample", "accession_number", "stage", "stageName", "sampleName")]

# timings of stages from the White et al. eLife (2017) publication of the data.
# time given in hours post-fertilization
timepoints <- data.frame(stage = unique(P_dswhite2017$stageName), 
                         hours_pf = c(0, .75, 2.25, 3, 4.3, 5.25, 6, 8, 10.3, 
                                      16, 19, 24, 30, 36, 48, 72, 96, 120),
                         stringsAsFactors = F, row.names = "stage")
P_dswhite2017$age <- timepoints[P_dswhite2017$stageName, "hours_pf"]
P_dswhite2017$batch <- factor(gsub(".*-(\\d)$", "\\1", P_dswhite2017$sampleName))

X_dswhite2017 <- X_dswhite2017[, P_dswhite2017$sample]

dswhite2017 <- list(g = X_dswhite2017, p = P_dswhite2017)
save(dswhite2017, file = paste0(data_folder, "dswhite2017.RData"), compress = "xz")

# cleanup
file.remove(g_file_dswhite2017)
rm(p_url_dswhite2017, g_url_dswhite2017, g_file_dswhite2017, X_dswhite2017, P_dswhite2017, timepoints)

`dslevin2016zeb`

geo_dslevin2016zeb <- "GSE60619"

g_url_dslevin2016zeb <- GEOquery::getGEOSuppFiles(geo_dslevin2016zeb, makeDirectory = FALSE, fetch_files = FALSE)
g_file_dslevin2016zeb <- paste0(data_folder, "dslevin2016zeb.txt.gz")
utils::download.file(url = as.character(g_url_dslevin2016zeb$url[2]), destfile = g_file_dslevin2016zeb)

X_dslevin2016zeb <- read.table(gzfile(g_file_dslevin2016zeb), h = T, sep = '\t', as.is = T, row.names = 1, comment.char = "")

# convert to tpm & ensembl_id
X_dslevin2016zeb <- X_dslevin2016zeb[rownames(X_dslevin2016zeb)%in%zeb_genes$ensembl_gene_id,]
X_dslevin2016zeb <- raw2tpm(rawcounts = X_dslevin2016zeb, 
                            genelengths = zeb_genes$transcript_length[match(rownames(X_dslevin2016zeb),
                                                                            zeb_genes$ensembl_gene_id)])


# pheno data
P_dslevin2016zeb <- Biobase::pData(GEOquery::getGEO(geo_dslevin2016zeb, getGPL = F)[[1]])

# filter relevant fields/samples
P_dslevin2016zeb <- P_dslevin2016zeb[, c("title", "geo_accession", "time (min after fertilization):ch1")]
colnames(P_dslevin2016zeb)[3] <- "time"
P_dslevin2016zeb$title <- as.character(P_dslevin2016zeb$title)

P_dslevin2016zeb <- P_dslevin2016zeb[P_dslevin2016zeb$title %in% colnames(X_dslevin2016zeb),]
X_dslevin2016zeb <- X_dslevin2016zeb[, P_dslevin2016zeb$title]

# formatting
P_dslevin2016zeb$title <- gsub('Metazome_ZF_timecourse_', '', P_dslevin2016zeb$title)
colnames(X_dslevin2016zeb) <- P_dslevin2016zeb$title

P_dslevin2016zeb$age <- as.numeric(P_dslevin2016zeb$time) / 60

dslevin2016zeb <- list(g = X_dslevin2016zeb, p = P_dslevin2016zeb)
save(dslevin2016zeb, file = paste0(data_folder, "dslevin2016zeb.RData"), compress = "xz")

# cleanup
file.remove(g_file_dslevin2016zeb)
rm(geo_dslevin2016zeb, g_url_dslevin2016zeb, g_file_dslevin2016zeb, X_dslevin2016zeb, P_dslevin2016zeb)

rm(zeb_genes, raw2tpm, fpkm2tpm)

Normalization & Quick look

dswhite2017$g <- limma::normalizeBetweenArrays(dswhite2017$g, method = "quantile")
dswhite2017$g <- log(dswhite2017$g + 1)

dslevin2016zeb$g <- limma::normalizeBetweenArrays(dslevin2016zeb$g, method = "quantile")
dslevin2016zeb$g <- log(dslevin2016zeb$g + 1)

dswhite2017$g[1:5, 1:5]
#>                    zmp_ph133_B zmp_ph133_D zmp_ph133_E zmp_ph133_F zmp_ph133_G
#> ENSDARG00000000001    2.192007    2.019082    1.929426    2.031762   1.9166338
#> ENSDARG00000000002    1.149510    1.188959    0.900076    1.185358   0.9783448
#> ENSDARG00000000018    2.456661    2.534134    2.224970    2.364784   2.5503750
#> ENSDARG00000000019    4.432509    4.529970    4.608232    4.533400   4.5923212
#> ENSDARG00000000068    4.406696    4.460862    4.267657    4.294028   4.1594844

head(dswhite2017$p, n = 5)
#>        sample accession_number         stage stageName sampleName age batch
#> 1 zmp_ph133_B       ERS1079239 Zygote:1-cell    1-cell   1-cell-1   0     1
#> 2 zmp_ph133_D       ERS1079240 Zygote:1-cell    1-cell   1-cell-2   0     2
#> 3 zmp_ph133_E       ERS1079241 Zygote:1-cell    1-cell   1-cell-3   0     3
#> 4 zmp_ph133_F       ERS1079243 Zygote:1-cell    1-cell   1-cell-4   0     4
#> 5 zmp_ph133_G       ERS1079244 Zygote:1-cell    1-cell   1-cell-5   0     5

Correlation Matrix

Plotting components

pca_dswhite2017 <- stats::prcomp(t(dswhite2017$g), rank = 25,
                                 center = TRUE, scale = FALSE)

Notice how the sampling is sparser towards the end of the time series, with dynamics being “wider”. Fitting splines on the components here will lead to a poor fit of the earlier timepoints.

To bypass this issue, we can use ranks instead of the timepoints.

# using data.table's rank function to get the "dense" tie method
dswhite2017$p$rank <- data.table::frank(dswhite2017$p$age, ties.method = "dense")

These dynamics will be fitted much more cleanly. To predict the data in a uniform time scale, we can just pick values on the rank scale such that they translate to a uniform series on the age scale with a simple linear warp, as will be done below.

Model fitting

Component number

nc <- sum(summary(pca_dswhite2017)$importance[3,] < .99) + 1
nc
#> [1] 67

Model

m_dswhite2017 <- ge_im(X = dswhite2017$g, p = dswhite2017$p, formula = "X ~ s(rank, bs = 'ds') + batch", nc = nc)

#>       aCC     aRE        MSE     aRMSE
#> 0.8593854 1.14166 0.03793716 0.1947746

Validation

Predict

We’ll be using a linear warp to get a uniform time series.

linwarp <- function(x, xyt, xc = 1, yc = 2){
  # Computes a linear interpolation of input x to y value
  # x = values of x to convert to y
  # xyt = table with known sets of x/y
  # xc, yc = column indices of x and y in xyt
  
  if(min(x) < min(xyt[,xc]) | max(x) > max(xyt[,xc]))
    stop("Some values of x are outside of the known x/y sets")
  
  # set up y to x conversion table
  xyt <- xyt[!duplicated.default(xyt[,xc]),]
  xyt <- xyt[order(xyt[,xc]),]
  
  xyt[,"dify"] <- c(0, diff(xyt[,yc]))
  xyt[,"difx"] <- c(1, diff(xyt[,xc]))
  xyt <- rbind(xyt[1,], xyt) # double 1st line for edge case

  xout <- unlist(sapply(x, function(xi){
    rsup <- which(xyt[-1,xc] >= xi)[1] + 1
    xyt[rsup-1, yc] + (xi - xyt[rsup-1, xc])/xyt[rsup, "difx"] * xyt[rsup, "dify"]
  }))
  
  return(xout)
}

# setup newdat
n.inter <- 200
newdat <- data.frame(age = seq(min(dswhite2017$p[, "age"]), max(dswhite2017$p[, "age"]), l = n.inter),
                     batch = rep("1", n.inter)) # predict as batch 1

# apply linwarp
newdat$rank <- linwarp(newdat$age, xyt = dswhite2017$p[, c("age", "rank")], xc = 1, yc = 2)

head(newdat)
#>         age batch     rank
#> 1 0.0000000     1 1.000000
#> 2 0.6030151     1 1.804020
#> 3 1.2060302     1 2.304020
#> 4 1.8090452     1 2.706030
#> 5 2.4120603     1 3.216080
#> 6 3.0150754     1 4.011596

# predict 
pred_m_dswhite2017 <- predict(m_dswhite2017, newdata = newdat)
pred_m_dswhite2017_comp <- predict(m_dswhite2017, newdata = newdat, as.c = TRUE)

We want a uniform series on the age scale, but have to input values on the rank scale in the model which is why we use linwarp(). To give a sense of what the function did, we can plot the ranks against the age.

Plot component predictions

On the rank scale :

Back on the age scale :

Build reference & stage samples

# make a 'reference object' 
r_dswhite2017 <- list(interpGE = pred_m_dswhite2017, time.series = newdat$age)

ae_dswhite2017 <- ae(dswhite2017$g, r_dswhite2017$interpGE, r_dswhite2017$time.series)
ae_dslevin2016zeb <- ae(dslevin2016zeb$g, r_dswhite2017$interpGE, r_dswhite2017$time.series)

On a model without using ranks

We’ll build the same model (not the optimal model !) without considering uneven sampling, for comparison. This will allow us to showcase some interpolation issues to look out for.

m_dswhite2017_norank <- ge_im(X = dswhite2017$g, p = dswhite2017$p, formula = "X ~ s(age, bs = 'ds') + batch", nc = nc)

#>       aCC       aRE      MSE    aRMSE
#> 0.8201534 0.6431004 2.668498 1.633554

pred_m_dswhite2017_norank <- predict(m_dswhite2017_norank, newdata = newdat)
pred_m_dswhite2017_comp_norank <- predict(m_dswhite2017_norank, newdata = newdat, as.c = TRUE)

Let’s plot the components on the rank and age scales, as before.

Back on the age scale :

We can already see that the model has trouble predicting the dynamics accurately. For example, we pick up noise in PC5 and flatten dynamics in PC6.

This has consequences on the estimates of the validation dataset, as you’ll see when we stage the samples.

# make a 'reference object' 
r_dswhite2017_norank <- list(interpGE = pred_m_dswhite2017_norank, time.series = newdat$age)

ae_dswhite2017_norank <- ae(dswhite2017$g, r_dswhite2017_norank$interpGE, r_dswhite2017_norank$time.series)
ae_dslevin2016zeb_norank <- ae(dslevin2016zeb$g, r_dswhite2017_norank$interpGE, r_dswhite2017_norank$time.series)

The “gaps” or “steps” you can see on the validation dataset’s estimates are due to interpolation bias : the picked up noise and flattened dynamics mentioned above. Essentially, the model errors create local “unlikely/unrealistic” gene expression zones in the interpolation, which will not correlate well with the samples of corresponding age. These zones will most often find themselves in between timepoints of the dataset used to build the reference, meaning the estimates of its samples are unaffected. However, the validation dataset has clear blank ranges of age estimates, around which are clustered the samples of corresponding age.

While we’ve used a suboptimal model here (which clearly displays model errors on component plots), some interpolation bias can be much more subtle. In such cases, this is something that can only be assessed using an external dataset.

References

Aeschimann, Florian, Pooja Kumari, Hrishikesh Bartake, Dimos Gaidatzis, Lan Xu, Rafal Ciosk, and Helge Großhans. 2017. “LIN41 Post-Transcriptionally Silences mRNAs by Two Distinct and Position-Dependent Mechanisms.” Molecular Cell 65 (3). Elsevier: 476–89.

Alter, Orly, Patrick O Brown, and David Botstein. 2000. “Singular Value Decomposition for Genome-Wide Expression Data Processing and Modeling.” Proceedings of the National Academy of Sciences 97 (18). National Acad Sciences: 10101–6.

Graveley, Brenton R, Angela N Brooks, Joseph W Carlson, Michael O Duff, Jane M Landolin, Li Yang, Carlo G Artieri, et al. 2011. “The Developmental Transcriptome of Drosophila Melanogaster.” Nature 471 (7339). Nature Publishing Group: 473.

Hendriks, Gert-Jan, Dimos Gaidatzis, Florian Aeschimann, and Helge Großhans. 2014. “Extensive Oscillatory Gene Expression During c. Elegans Larval Development.” Molecular Cell 53 (3). Elsevier: 380–92.

Levin, Michal, Leon Anavy, Alison G Cole, Eitan Winter, Natalia Mostov, Sally Khair, Naftalie Senderovich, et al. 2016. “The Mid-Developmental Transition and the Evolution of Animal Body Plans.” Nature 531 (7596). Nature Publishing Group: 637.

White, Richard J, John E Collins, Ian M Sealy, Neha Wali, Christopher M Dooley, Zsofia Digby, Derek L Stemple, et al. 2017. “A High-Resolution mRNA Expression Time Course of Embryonic Development in Zebrafish.” Elife 6. eLife Sciences Publications Limited: e30860.

RAPToR - Building References

Romain Bulteau

mai 2020

Preamble

The data

Selecting / Preparing a dataset

Databases

What to look out for

ID formatting

Normalize and log expression

Observing the data

Correlation

Plots

Code

Plotting components

Plot

Code

Plotting random genes

Plots

Code

The gene expression interpolation model

About the model

The GEIM interface

Finding the appropriate model and parameters

Model type

Model performance

Plots

Code

Number of components

Comparing formulas

Predicting new timepoints

Validating / Checking model predictions

Checking model predictions against components

Plots

Code

Plots

Code

Staging samples

Plots

Code

Reference-Building examples

C. elegans larval development

The data

Hide

Show

dsaeschimann2017

dshendriks2014

Normalization & Quick look

Correlation Matrix

Plotting components

Model fitting

Component number

Model

Validation

Predict

Plot component predictions

Build reference & stage samples

D. melanogaster embryonic development

The data

Hide

Show

dsgraveley2011

dslevin2016dmel

Normalization & Quick look

Correlation Matrix

Plotting components

Model fitting

Component number

Model

Validation

Predict

Plot component predictions

Build reference & stage samples

Danio rerio embryonic development

The data

Hide

Show

dswhite2017

dslevin2016zeb

Normalization & Quick look

`RAPToR` - Building References

`dsaeschimann2017`

`dshendriks2014`

`dsgraveley2011`

`dslevin2016dmel`

`dswhite2017`

`dslevin2016zeb`